/* Stata code for Regional Heterogeneity in Environmental Quality: The Role of Firm Production 
   Networks and Trade (Park, Howard, and Ridley; Journal of the Association of Environmental
   and Resource Economists) */

/* This do file combines the various data sources on states' emissions/emissions intensities,
   market access measures (HRI), and intermediate trade to generate Figures 1 through 3a in
   the stylized facts section.
   
   See the manuscript for detailed descriptions of variables and their sources.
   
   This do file requires the Stata package grc1leg */

*net install grc1leg.pkg

/* Indicate the directory where all of the data is stored and figures will be outputed */

local directory = "C:\users\wridley\desktop\Stata"

cd "`directory'"


***************************************************************************
****** Figure 1 - Rank of states by number of establishments and HRI ******
***************************************************************************

use "establishments_and_HRI_by_state.dta", clear

gen byte core25 = mean_HRI > 1.33 // Indicator for states in the top 25th percentile of HRI (i.e., core states)

/* Make scatter plots of state-level rankings in establishment counts and HRI index */
egen rank_estbcount = rank(estb), field
egen rank_hriindex = rank(mean_HRI), field

corr rank_hriindex rank_estbcount // 0.8111
local corr = string(r(C)[2,1],"%03.2f")

graph set eps fontface "Arial"

tw (function y = x, ra(rank_estbcount) clpat(dash)) ///
	(sc rank_estbcount rank_hriindex if core25 == 1, ///
	mc(red) mlabel(state) msymbol(s) mlabcolor(black) mlabsize(vsmall))  ///
	(sc rank_estbcount rank_hriindex if core25 == 0, ///
	mc(blue) mlabel(state) msymbol(Oh) mlabcolor(black) mlabsize(vsmall)), ///
	title("") xsize(6) ysize(6) ///
	xscale(reverse) yscale(reverse) ///
	xmtick(#48) ymtick(#48) ///
	ytitle("Rank, Number of Establishments") xtitle("Rank, HRI") ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "45{c 176} line") pos(6)) ///
	ylabel(1 12 24 36 48,angle(h) nogrid) ///
	xlabel(1 12 24 36 48, nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	text(45 5 "{it:R} = `corr'", size(small) color(gs10))
	
*graph export "fig01_establishments_HRI_scatter.eps", as(eps) replace


***********************************************************************
****** Figure 2 - Emissions/emissions intensity and HRI by state ******
***********************************************************************

use "emissions_and_HRI_by_state.dta", clear

/* Calculate emission intensity by state variable (kiloton of local pollutant / billion chained 2009 dollars) */
gen eint = emissions / output

gen byte core25 = mean_HRI > 1.33 // Indicator for states in the top 25th percentile of HRI (i.e., core states)

/* Generating log emissions and emissions intensities and calculating correlations with HRI */
foreach i in CO PM25 NOx {

	gen ln_`i' = log(emissions) if type == "`i'"			
	gen ln_`i'eint = log(eint) if type == "`i'"

	di "`i'"
	corr emissions mean_HRI if type == "`i'"
	corr eint mean_HRI if type == "`i'"

	corr ln_`i' mean_HRI if type == "`i'"
	local corr_`i' = string(r(C)[2,1],"%03.2f")
	
	corr ln_`i'eint mean_HRI if type == "`i'"
	local corr_eint_`i' = string(r(C)[2,1],"%03.2f")
	
}

/* Manually defining positions for marker labels */

qui foreach ii in a b c d e f {
	
	capture gen mlabelpos_`ii' = .
	replace mlabelpos_`ii' = 3 if state == "CA"
	replace mlabelpos_`ii' = 3 if state == "NJ"
	replace mlabelpos_`ii' = 3 if state == "NY"
	
}

qui {

	replace mlabelpos_a = 8 if state == "FL"
	replace mlabelpos_b = 8 if state == "FL"
	replace mlabelpos_c = 7 if state == "FL"
	replace mlabelpos_d = 6 if state == "FL"
	replace mlabelpos_e = 6 if state == "FL"
	replace mlabelpos_f = 12 if state == "FL"
	
	replace mlabelpos_f = 2 if state == "GA"

	replace mlabelpos_a = 9 if state == "IL"
	replace mlabelpos_b = 6 if state == "IL"
	replace mlabelpos_c = 3 if state == "IL"
	replace mlabelpos_d = 4 if state == "IL"
	replace mlabelpos_f = 6 if state == "IL"

	replace mlabelpos_a = 1 if state == "KY"
	replace mlabelpos_b = 9 if state == "KY"
	replace mlabelpos_c = 3 if state == "KY"
	replace mlabelpos_d = 2 if state == "KY"
	replace mlabelpos_f = 1 if state == "KY"

	replace mlabelpos_a = 2 if state == "MI"
	replace mlabelpos_d = 2 if state == "MI"

	replace mlabelpos_a = 9 if state == "SC"
	replace mlabelpos_b = 6 if state == "SC"
	replace mlabelpos_c = 12 if state == "SC"
	replace mlabelpos_d = 12 if state == "SC"
	replace mlabelpos_e = 9 if state == "SC"
	replace mlabelpos_f = 1 if state == "SC"

	replace mlabelpos_a = 9 if state == "TN"
	replace mlabelpos_b = 12 if state == "TN"
	replace mlabelpos_c = 7 if state == "TN"
	replace mlabelpos_d = 9 if state == "TN"
	replace mlabelpos_e = 12 if state == "TN"
	replace mlabelpos_f = 7 if state == "TN"

	replace mlabelpos_a = 2 if state == "TX"
	replace mlabelpos_b = 12 if state == "TX"

	replace mlabelpos_a = 3 if state == "WA"
	replace mlabelpos_c = 3 if state == "WA"
	replace mlabelpos_d = 2 if state == "WA"
	replace mlabelpos_e = 4 if state == "WA"

}

graph set eps fontface "Arial"

/* Individual scatter plots of total emissions by pollutant/state and HRI */

qui {
/* (a) CO Emissions */	
reg ln_CO mean_HRI	// linear fit of log(CO) and HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)

tw  (lfit ln_CO mean_HRI if type == "CO") ///
	(sc ln_CO mean_HRI if type == "CO" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_a)) ///
	(sc ln_CO mean_HRI if type == "CO" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(a) CO Emissions", color(black) size(medium)) ///
	ytitle("log(Total Emission)") xtitle("HRI") ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-3 3.7 "y = `beta'*HRI + `cons'" "{it:R} = `corr_CO'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logCOemission_HRI_2017, replace) ///
	nodraw
	
/* (b) PM2.5 Emissions */
reg ln_PM25 mean_HRI	// linear fit of log(PM2.5) and HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)
	
tw  (lfit ln_PM25 mean_HRI) ///
	(sc ln_PM25 mean_HRI if type == "PM25" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_b)) ///
	(sc ln_PM25 mean_HRI if type == "PM25" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(b) PM{sub:2.5} Emissions", color(black) size(medium)) ///
	ytitle("log(Total Emission)") xtitle("HRI") ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-3.25 3.7 "y = `beta'*HRI `cons'" "{it:R} = `corr_PM25'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logPM25emission_HRI_2017, replace) ///
	nodraw
	
/* (c) NOx Emissions */
reg ln_NOx mean_HRI	// linear fit of log(NOx) and HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)
	
tw  (lfit ln_NOx mean_HRI) ///
	(sc ln_NOx mean_HRI if type == "NOx" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_c)) ///
	(sc ln_NOx mean_HRI if type == "NOx" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(c) NO{sub:x} Emissions", color(black) size(medium)) ///
	ytitle("log(Total Emission)") xtitle("HRI") ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-4 3.7 "y = `beta'*HRI + `cons'" "{it:R} = `corr_NOx'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logNOxemission_HRI_2017, replace) ///
	nodraw

/* Individual scatter plots of emission intensities by pollutant/state and HRI */

/* (d) CO Emission Intensity */
reg ln_COeint mean_HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)
	
tw  (lfit ln_COeint mean_HRI) ///
	(sc ln_COeint mean_HRI if type == "CO" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_d)) ///
	(sc ln_COeint mean_HRI if type == "CO" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(d) CO Emission Intensity", color(black) size(medium)) /// 
	ytitle("log(Emission Intensity)") xtitle("HRI") ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-7 3.7 "y = `beta'*HRI `cons'" "{it:R} = `corr_eint_CO'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logCOintensity_HRI_2017, replace) ///
	nodraw

/* (e) PM2.5 Emission Intensity */
reg ln_PM25eint mean_HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)
	
tw  (lfit ln_PM25eint mean_HRI) ///
	(sc ln_PM25eint mean_HRI if type == "PM25" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_e)) ///
	(sc ln_PM25eint mean_HRI if type == "PM25" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(e) PM{sub:2.5} Emission Intensity", color(black) size(medium)) ///
	ytitle("log(Emission Intensity)") xtitle("HRI") ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-7.25 3.7 "y = `beta'*HRI `cons'" "{it:R} = `corr_eint_PM25'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logPM25intensity_HRI_2017, replace) ///
	nodraw 
	
/* (f) NOx Emission Intensity */
reg ln_NOxeint mean_HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)
	
tw  (lfit ln_NOxeint mean_HRI) ///
	(sc ln_NOxeint mean_HRI if type == "NOx" & core25 == 1, ///
	mc(red) mlabel(state) mlabc(black) msymbol(s) mlabsize(vsmall) mlabgap(0.3) mlabvpos(mlabelpos_f)) ///
	(sc ln_NOxeint mean_HRI if type == "NOx" & core25 == 0, ///
	mc(blue) mlabel("") msymbol(Oh)), ///
	title("(f) NO{sub:x} Emission Intensity", color(black) size(medium)) ///
	ytitle("log(Emission Intensity)") xtitle("HRI") ///
	ylabel(0(2)-10,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	legend(region(lcolor(white)) rows(1) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line") size(small)) ///
	text(-9 3.7 "y = `beta'*HRI `cons'" "{it:R} = `corr_eint_NOx'", size(small) color(gs10)) ///
	xsize(6) ysize(6) ///
	saving(logNOxintensity_HRI_2017, replace) ///
	nodraw
	
}
	
grc1leg "logCOemission_HRI_2017.gph" "logPM25emission_HRI_2017.gph" "logNOxemission_HRI_2017.gph" ///
	"logCOintensity_HRI_2017.gph" "logPM25intensity_HRI_2017.gph" "logNOxintensity_HRI_2017", ///
	leg("logCOemission_HRI_2017.gph") ring(6) rows(2) ///
	imargin(small) ///
	graphr(style(none) color(white))
	
*graph export "fig02_emissions_HRI_scatter.eps", as(eps) replace
	
/* Cleaning up individual graph files */
rm logCOemission_HRI_2017.gph
rm logCOintensity_HRI_2017.gph
rm logPM25emission_HRI_2017.gph
rm logPM25intensity_HRI_2017.gph
rm logNOxemission_HRI_2017.gph
rm logNOxintensity_HRI_2017.gph


******************************************************************
****** Figure 3 - Net intermediate imports and HRI by state ******
******************************************************************

use "intermediate_trade_and_HRI_by_state.dta", clear

gen byte core25 = mean_HRI > 1.33 // top 25 percentile
gen netimports = imports - exports 	// net intermediate imports by state

corr netimports mean_HRI	// calculating correlation between net intermediate imports and HRI
local corr = string(r(C)[2,1],"%03.2f")

drop if state == "CA" // drop visual outliers for better illustration; including CA doesn't change the relationship

reg netimports mean_HRI	// linear fit of netimports and HRI
	local beta = string(_b[mean_HRI],"%03.2f")
	local cons = subinstr(string(_b[_cons],"%03.2f"),"-","- ",.)

graph set eps fontface "Arial"

tw (lfit netimports mean_HRI) /// 
	(sc netimports mean_HRI if core25 == 1, ///
	mc(red) mlabel(state) msymbol(s) mlabcolor(black)) ///
	(sc netimports mean_HRI if core25 == 0, /// 
	mc(blue) mlabel("") msymbol(Oh) mlabcolor(black)), ///
	title("") ytitle("Net Intermediate Imports (billion USD)") xtitle("HRI") ///
	yscale(titlegap(5)) ///
	ylabel(,angle(h) nogrid) ///
	xlabel(,nogrid) ///
	graphr(style(none) color(white)) ///
    plotr(lc(black) lw(medthin) m(l=2 r=2)) ///
	text(-7.75 3.1 "y = `beta'*HRI `cons'" "{it:R} = `corr'", size(medsmall) color(gs10)) ///
	xsize(6) ysize(6) ///
	legend(region(lcolor(white)) pos(6) rows(1) colgap(15) order(2 "Core States" 3 "Periphery States" 1 "Fitted Line"))
	
*graph export "fig03a_scatterplot_USAtrade_netimports_HRI.eps", as(eps) replace